function [M,var_M,Mhatb] = compute_data_moments_IRF_cf_persist(X_monthly,Nu)
 
% Periods 1-6: pre shock
% Periods 7-14: post shock

% Replications of the bootstrap
B                       = Nu.B;

% Point estimate
X                       = X_monthly(7:end,:);
TimeZeroX               = repmat(X_monthly(1,:),[8,1]);
dX                      = X - TimeZeroX;
M                       = mean(dX,2);

% Bootstrapped covariance matrix
Mb                      = NaN(numel(M),B);
for b=1:B
   dXb                      = dX(:,Nu.indices1(:,b)); 
   Mb(:,b)                  = mean(dXb,2);
end
Mhatb                   = mean(Mb,2);
Dummy                   = Mb-repmat(Mhatb,[1 B]);
var_M                   = diag(diag(1/(B-1)*Dummy*(Dummy')));

% Diagonal covariance matrix
var_M                   = (1e-5)*eye(numel(M));

end

